Cite this Article
Lu Jia, Zhou Huaichun. Finite-difference time-domain modeling of curved material interfaces by using boundary condition equations method. Chinese Physics B, 2016, 25(9): 090203  
Permissions
Finite-difference time-domain modeling of curved material interfaces by using boundary condition equations method
Lu Jia, Zhou Huaichun†,
Department of Thermal Engineering, Tsinghua University, Beijing 100084, China

 

† Corresponding author. E-mail: hczh@mail.tsinghua.edu.cn

Project supported by the National Natural Science Foundation of China (Grant No. 51025622).

Abstract
Abstract

To deal with the staircase approximation problem in the standard finite-difference time-domain (FDTD) simulation, the two-dimensional boundary condition equations (BCE) method is proposed in this paper. In the BCE method, the standard FDTD algorithm can be used as usual, and the curved surface is treated by adding the boundary condition equations. Thus, while maintaining the simplicity and computational efficiency of the standard FDTD algorithm, the BCE method can solve the staircase approximation problem. The BCE method is validated by analyzing near field and far field scattering properties of the PEC and dielectric cylinders. The results show that the BCE method can maintain a second-order accuracy by eliminating the staircase approximation errors. Moreover, the results of the BCE method show good accuracy for cylinder scattering cases with different permittivities.

PACS: 02.60.–x;02.70.–c;41.20.Jb;46.15.–x
Keyword:finite-difference time-domain;curved surface;staircase error;boundary condition equation
1. Introduction

The interactions of electromagnetic waves with a curved surface have been widely discussed in many applications, such as material science, optics, photonic crystals, antenna design, and remote sensing.[1–6] In recent decades, a number of numerical methods have been reported for the interactions of electromagnetic waves with curved surfaces, including finite difference methods,[7,8] finite element methods,[9,10] and integral equation methods.[11,12] Among these methods, the finite-difference time-domain method is widely used in computational electrodynamics for easy implementation and accurate results.[13–17] When modeling a curved surface, however, the numerical accuracy may be reduced, because the standard FDTD algorithm uses staircase grids to approximate the material interface.[18] The staircase approximation, as depicted in Fig. 1, is generally adopted because of the discretization of the surface position. Obviously, the scattering properties of the curved surface will change going from the real shape to the staircase shape, especially if the mesh is coarse and the dielectric contrast is large.[19]

Fig. 1. Staircase approximation of a circular cylinder in FDTD calculation. Here, UPML means uniaxial perfectly matched layer. TF/SF: total field/scattering field.

A possible solution to the staircase approximation problem is using non-orthogonal grids or curvilinear coordinates to approximate the curved surface.[20,21] Obviously, while improving the numerical accuracy, these methods would increase the complexity of the algorithm and cause numerical artifacts.[19] A more commonly used method is the conformal FDTD method.[22] However, the numerical stability of the contour-path FDTD method may not be guaranteed, because there are electric nodes whose values are borrowed from neighbor cells.[23] Research is ongoing for the improvements of the contour-path FDTD method based on the numerical instability problems.[23,24]

Another possible way for solving the staircase approximation problem is using the effective permittivity to model the dielectric interface. There are several works presenting different kinds of effective permittivity methods, including the volume average method, a first-neighbor average method, the Maxwell–Garnett method, an inverted Maxwell–Garnett method, and the Bruggeman formula method.[19] The problem of the effective permittivity method resides in the selection of the most suitable effective permittivity value. There have been improvements to the effective permittivity method.[25,26] Among these improvements, the sub-pixel smoothing method shows second-order accuracy by employing the mean permittivity for surface-parallel electric components and the harmonic mean permittivity for surface-perpendicular electric components.[27,28]

A third approach is to simulate the real shape by considering the physical location of curved surfaces.[29,30] These methods show good accuracy and almost need no additional computational cost. However, the complexity of the algorithm would be increased due to the modification of the Yee algorithm. According to the classical electromagnetic theory, material interface boundary condition equations are essential for the scattering problems calculation.[31] Therefore, we utilize the Yee algorithm as usual and treat the surface interface as a kind of interface boundary. The layer of grids across the boundary is calculated by considering the boundary condition equations and physical location of the surface interface. Thus, while maintaining the simplicity and computational efficiency of the Yee scheme, the algorithm can solve the staircase approximation problem.

In this work, the boundary condition equations (BCE) method is utilized to deal with the staircase approximation problem in curved surface scattering. The paper is organized as follows. Firstly, the two-dimensional BCE algorithm is proposed in this paper. Subsequently, the BCE method is validated by analyzing the induced surface current on a PEC circular cylinder. Furthermore, the near field scattering properties of the dielectric circular cylinders are investigated, and the relative errors of different methods are compared. Then, we discuss the far field scattering properties of dielectric cylinders and show the effect of different dielectric permittivities on different methods. Finally, some concluding remarks are derived.

2. Numerical method

In this section, the BCE method is introduced to solve the staircase approximation problem in curved surface scattering.

2.1. Two-dimensional FDTD method

The finite difference time domain method is a state-of-the-art method for solving Maxwell's equations in complex geometries.[13] Considering linear, isotropic, non-dispersive, and lossy media, the FDTD method solves Maxwell’s curl equations as follows:[13]

where E and H are electric and magnetic field vectors, Jsource and Msource are electric and magnetic sources, ε, μ, σ, and σ* are electric permittivity, magnetic permeability, electric conductivity, and equivalent magnetic conductivity.

For TE wave incidence, Maxwell’s equations in two-dimensional FDTD can be found in previous publications as[13]

where

Two-dimensional FDTD equations for TM wave incidence can also be found in previous publications.[13]

2.2. Dielectric interface boundary condition

The FDTD computation setup is presented in Fig. 1. The FDTD calculation region is surrounded by the uniaxial perfectly matched layer (UPML) to simulate electromagnetic wave propagation outside the computation region without reflection. In the computation region, there are plane waves incident on the infinite cylinder. The plane waves are added by using the total field/scattering field (TF/SF) method.[13] As shown in Fig. 1, the FDTD computation region is divided by discrete grids. Within the grids, there are Hz, Ex, and Ey components represented by circles and arrows. There are two regions painted by different colors. The two regions represent two kinds of materials. The red line between the two materials is the material interface, which is usually approximated by the black line with a staircase shape in FDTD simulations. By using the staircase approximation, electromagnetic wave scattering will change going from the real shape to the staircase shape, especially if the mesh is coarse and the dielectric contrast is large.[19] In order to make the material parameters change smoothly and improve the numerical accuracy, the averaged material parameters of adjacent grids have been used by earlier researchers. Usually, there are two kinds of method, namely arithmetic mean and harmonic mean method,[32] they are

Arithmetic mean:[32]

Harmonic mean:[32]

However, these methods still suffer from the drawbacks of staircase approximation. To improve the numerical accuracy, the boundary condition equations on the material interface must be considered. For the curved surface case, we can use the same Yee algorithm as usual and just apply the boundary condition equations on the interface boundary. To calculate the grids across the interface boundary, the physical location of the interface boundary must be considered rigorously.[29]

Figure 2 shows the electromagnetic fields configuration near the material interface. Consider the scheme for updating the Hz field, the magnetic field point and the nearest four electric field points are treated as a group. When an interface passes through the group of points, it separates some electric fields and the magnetic field into different media. In Fig. 2, for example, the material interface intersects one line connecting Ey(i, j + 1/2) and Hz(i + 1/2, j + 1/2) at point Int1 and another line connecting Ex(i + 1/2, j + 1) and Hz(i + 1/2, j + 1/2) at point Int2. Thus, Ex(i + 1/2, j + 1) and Ey(i, j + 1/2) are separated from the other three components. Here, φx denotes the ratio of the distance between point Int1 and Hz(i + 1/2, j + 1/2) to the distance between Hz(i – 1/2, j + 1/2) and Hz(i + 1/2, j + 1/2). φy denotes the ratio of the distance between point Int2 and Hz(i + 1/2, j + 1/2) to the distance between Hz(i + 1/2, j + 3/2) and Hz(i + 1/2, j + 1/2). ET(int11) and EN(int11) are the tangential and normal components of the electric field respectively at point Int1; ET(int21) and EN(int21) are those components at point Int2. In addition, it can be seen from Fig. 2 that the surface interface makes an angle of α with the x direction.

Fig. 2. Electromagnetic fields configuration for the material interface. Here, ET and EN are the tangential and normal components of the electric field.

To update Hz(i + 1/2, j + 1/2) in Fig. 2, both and in Eq. (2c) must be replaced by new values of and in medium 2. The new in medium 2 can be approximated as

i.e.,

Here, int12 means that the point is just below the boundary interface in media 2. For grids across the boundary, coefficients of Ca, Cb, Da, and Db in Eq. (2) are calculated based on the location of the boundary. To determine the location of the boundary, the surface location, intersection ratio, and surface slope angle are calculated in the preprocessing stage.

Similarly, the formula of the new is

where and in Eqs. (4)–(6) are unknown parameters, which can be calculated by the boundary condition equations near the material interface. Boundary condition equations at the dielectric interface state that the tangential component of the electric field is continuous across the boundary, whereas the normal component of the electric field changes discontinuously across the boundary.[31] The calculation of will be shown by applying boundary condition equations.

According to the forward difference approximation with the first-order accuracy, we obtain

and based on the forward difference approximation with the second-order accuracy, we obtain

Substituting Eq. (8) into Eq. (7), we can obtain

and then, can be approximated as

i.e.,

Here, int21 means that the point is just above the boundary interface in media 1.

Similarly, the formula of is

With and known, we can derive and using coordinates transformation

To calculate and , the boundary condition equations at the dielectric interface should be applied as[33]

where D = εE and B = μH are the electric flux density and magnetic flux density.

Then, and can be calculated as

In addition,

Therefore, the boundary condition equations for material interface can be simplified as

where the matrices R and T are defined as

Thus, we derive the formulas of based on boundary condition equations. Similarly, can be calculated. Then, the Hz(i + 1/2, j + 1/2) field is updated according to Eq. (2c). Here, new values of and in medium 2 are utilized in Eq. (2c). Parameters of medium 2 are used in Da and Db.

Consider the scheme of updating the Ex(i + 1/2, j + 1) field, the electric field point and the nearest two magnetic field points are treated as a group. When an interface passes through the group of points, it separates one magnetic field and the electric field into different media. The configuration can be presented in Fig. 2.

To update Ex(i + 1/2, j + 1) in Fig. 2, in Eq. (2a) must be replaced by a new value of in medium 1.

i.e.,

Here, is an unknown parameter, which can be calculated by applying the boundary condition equations in Eq. (14) as

i.e.,

The scheme for updating the Ey(i, j + 1/2) field in Fig. 2 is the same as that of the Ex(i + 1/2, j + 1) field, in Eq. (2b) is replaced by the new value of in medium 1.

2.3. PEC interface boundary condition

In the case of a PEC interface, the boundary condition equations are different from that of the dielectric interface. Assuming media 2 in Fig. 2 is PEC, the boundary condition equations become[33]

Then, Hz(i + 1/2, j + 1/2) = 0.

To update Ex(i + 1/2, j + 1) in Fig. 2, in Eq. (2a) must be replaced by a new value of in medium 1.

i.e.,

Here, is an unknown parameter, which can be calculated by applying the boundary condition equations in Eq. (23),

i.e.,

The scheme for updating the Ey(i, j + 1/2) field in Fig. 2 is the same as that of the Ex(i + 1/2, j + 1) field, in Eq. (2b) is replaced by the new value of in medium 1.

Thus, we have set up the methodology for adding the boundary condition equations in a two-dimensional TE FDTD algorithm. In summary, the flow chart of the BCE method is drawn in Fig. 3. Similarly, we can derive the methodology for adding the boundary condition equations in the two-dimensional TM FDTD algorithm.

Fig. 3. Flow chart of the boundary condition equations method.

The BCE method utilizes the standard FDTD algorithm with a second-order accuracy. All the approximations across the boundary are based on the one side difference approximation or non-centered difference approximation. Thus, it is sufficient to guarantee a second-order accuracy globally.[30] In addition, the BCE method can be extended to solve three-dimensional problems. The surface location (i, j, k) and intersection ratio (φx, φy, φz) should be calculated instead of the surface location (i, j) and intersection ratio (φx, φy) in the preprocessing stage, also, matrix R and T must be replaced by a coordinate transform matrix in three dimensions. The boundary condition equations are similar to those of two-dimensional problems.

3. Results and discussion

In this section, analytical and numerical data validations of the BCE method are presented. Comparisons of different methods are made, including the staircase approximation (SA) method, arithmetic mean (AM) method, analytical solution (AS) method, and BCE method. Results from the sub-pixel smoothing (SPS) method are also given for method validation.[27,28]

3.1. PEC circular cylinder

To validate the FDTD algorithm, the induced surface current of a PEC circular cylinder is analyzed. The PEC circular cylinder has a radius of a = 5/ki and is under transverse electric (TE) illumination. Here, ki = 2π/λ is the wave vector. Figure 4 is a plot of the BCE method predicted surface currents and the SA method predicted surface currents compared to the reference data.[22] It can be seen from Fig. 4 that the BCE method predicted surface currents agree well with the series solution and contour FDTD results, whereas the SA method predicted surface currents do not agree well with the reference data.

Fig. 4. Surface induced currents comparison. Series solution and contour FDTD are reference data. SA and BCE are staircase approximation and boundary condition equations method.
3.2. Dielectric circular cylinder, near field

In this part, the near field scattering properties of a kia = 5 dielectric circular cylinder is analyzed. The permittivity and permeability in free space are given as ε0 = 1, μ0 = 1, σ0 = 0, and σm,0 = 0. The material parameters of the dielectric circular cylinder are εd = 2, μd = 1, σd = 0, and σm,d = 0. The incident waves are time-harmonic TE mode fields with an amplitude of |Hin| = 1.0 A/m and transverse magnetic (TM) mode fields with an amplitude of |Ein| = 1.0 V/m. The incident waves have a wavelength of λ = 1 m. The grid size is chosen as λ/160, and the simulation time is chosen to be t = 1.67 × 10−7 s to make sure that the numerical results are convergent. We use the SA, BCE, and AS methods to calculate a total field region of [–2 : 2] m×[–2 : 2] m. For the infinite cylinder scattering cases, analytical solutions can be derived from previous publications.[33–35]

To compare the local error in the whole computation region, the normalized L2 error for TE wave incidence can be defined as

The normalized L2 error for TM wave incidence is similar to that for TE wave incidence.

Figure 5 shows normalized L2 errors of the BCE method predicted total field amplitude and the SA and AM methods-predicted total field amplitude. It can be seen from Fig. 5 that the SA, AM, and BCE methods are convergent with the decrease of the grid size for both TE and TM illuminations. The BCE method can maintain a second-order accuracy by eliminating the staircase approximation error, whereas the SA and AM method can only achieve a first-order accuracy. For the TE wave incidence, the BCE method with a grid size of λ/40 is good enough to achieve a better accuracy than the SA method with a grid size of λ/160. For the TM wave incidence, the SA and AM methods have the same normalized L2 error, because there are no magnetic permeability differences between the circular cylinder and air.

Fig. 5. The normalized L2 errors for plane wave incidence on the dielectric cylinder. SA, AM, and BCE are staircase approximation, arithmetic mean, and boundary condition equations method, respectively. (a) TE wave illumination and (b) TM wave illumination.
3.3. Multilayer dielectric cylinder, near field

In addition, the BCE method can be used to deal with multilayer object scattering problems. For the multilayer circular cylinder scattering cases, analytical solutions can be derived based on the Mie theory.[36,37] In this part, the near field scattering properties of two double-layer dielectric circular cylinders is analyzed. The radius of the first double-layer circular cylinder is a = 5/(2π) m and ai = 5/(4π) m. The radius of the second double-layer circular cylinder are a = 15/(2π) m and ai = 15/(4π) m. The material parameters of the outer dielectric circular cylinders are εd = 2, μd = 1, σd = 0, and σm,d = 0. The material parameters of the inner dielectric circular cylinders are εd,i = 4, μd = 1, σd,i = 0, and σm,d,i = 0. The incident waves have a wavelength of λ = 1 m. We use the SA, BCE, SPS, and AS methods to calculate a total field region of [–2 : 2] m×[–2 : 2] m.

Fig. 6. The normalized L2 errors for plane wave incidence on the double layer dielectric cylinder (a = 5/(2π), ai = 5/(4π)). SA, BCE, and SPS are staircase approximation, boundary condition equations, and sub-pixel smoothing method, respectively. (a) TE wave illumination and (b) TM wave illumination.

Figures 6 and 7 show normalized L2 errors of the BCE method predicted total field amplitude and the SA and SPS methods-predicted total field amplitude. It can be seen from the two figures that the BCE and SPS methods can maintain a second-order accuracy by eliminating the staircase approximation error, whereas the SA method can only achieve a first-order accuracy. Meanwhile, it can be seen from the two figures that the results of the BCE method show better accuracy than that of the SPS method. The reason is that the BCE method is able to consider the physical location of curved surfaces according to the surface location and intersection ratio. Comparing Fig. 6 with Fig. 7, it can be seen that the results of the BCE and SPS methods show less improvements of the staircase errors for larger structures. The reason is that the dielectric interface becomes smoother with the increase of the structure size.

Fig. 7. The normalized L2 errors for plane wave incidence on the double layer dielectric cylinder (a = 15/(2π), ai = 15/(4π)). SA, BCE, and SPS are staircase approximation, boundary condition equations, and sub-pixel smoothing methods, respectively. (a) TE wave illumination and (b) TM wave illumination.
3.4. Dielectric circular cylinder, far field

To compare the performance of the SA method and BCE method, the far field RCS of a dielectric circular cylinder with a = 0.2λ is analyzed. The permittivity and permeability in free space are given as ε0 = 1, μ0 = 1, σ0 = 0, and σm,0 = 0. The material parameters of the dielectric circular cylinder are εd = 9, μd = 1, σd = 0, and σm,d = 0. The incident waves are time-harmonic TE mode fields with an amplitude of |Hin| = 1.0 A/m and TM mode fields with an amplitude of |Ein| = 1.0 V/m. The incident waves have a wavelength of λ = 1 m. The grid size is chosen as λ/80, and the simulation time is chosen to be t = 1.25 × 10−8 s to make sure that the numerical results are convergent. We use the SA, BCE, and analytical solution methods to calculate the far field RCS.

Figure 8 shows the BCE method predicted RCS and the SA method predicted RCS compared with the analytical solution results. Obviously, the results of the BCE method agree well with those of the analytical solution method for large dielectric contrast cases. It can be seen from Fig. 8 that the results of the SA method with a grid size of λ/80 do not agree well with the analytical solution results for the large dielectric contrast case. The reason is that the scattering properties of a curved surface will change going from the real shape to the staircase shape, especially if the mesh is coarse and the dielectric contrast is large.[19] Meanwhile, the results of the SA method show better agreement with the analytical solution results with a smaller grid size of λ/320.

Fig. 8. Far field RCS comparison (a = 0.2λ, εd = 9). SA and BCE are staircase approximation and boundary condition equations methods. (a) TE wave illumination and (b) TM wave illumination.

To demonstrate the method capability, the far field RCS of two cylinders with complex geometries are analyzed. The two cylinders are presented in Figs. 9 and 10. The cylinder in Fig. 9 has a square shape with a = λ, and the cylinder in Fig. 10 is a combination of a circular cylinder and a rectangle cylinder. The material parameters of the two cylinders are εd = 9, μd = 1, σd = 0, and σm,d = 0. The incident waves have a wavelength of λ = 1 m. The grid size is chosen as λ/80, and the simulation time is chosen to be t = 1.67 × 10−7 s to make sure that the numerical results are convergent. We use the SA and BCE methods to calculate the far field RCS, and results of the MOM method are used as the benchmark for the comparison with the FDTD results.[35]

Fig. 9. Far field RCS comparison (a = λ, εd = 9). SA and BCE are staircase approximation and boundary condition equations methods. (a) TE wave illumination and (b) TM wave illumination.
Fig. 10. Far field RCS comparison (a = 0.2λ, εd = 9). SA and BCE are staircase approximation and boundary condition equations methods. (a) TE wave illumination and (b) TM wave illumination.

Figures 9 and 10 show the BCE method predicted RCS and the SA method predicted RCS compared with the results of the MOM method. Obviously, the results of the BCE method agree well with those of the MOM method for both TE and TM wave illuminations, whereas the results of the SA method with a grid size of λ/80 do not agree well with that of the MOM method for TE wave illumination. The reason is that the sloped surfaces on the square cylinder will introduce a staircase error in the FDTD simulation.[19] The BCE method is able to eliminate the staircase approximation error and restore the second-order accuracy in FDTD simulations.

4. Conclusion

The BCE method is proposed to deal with the curved surface scattering in the FDTD calculation. In this method, the curved surface interfaces are modeled by adding the boundary condition equations. Thus, the real shapes of a curved surface can be simulated and the staircase approximation errors can be eliminated. To validate the BCE method, various scattering properties are analyzed, including the induced surface current of a PEC circular, and the near field and far field scattering of the dielectric cylinders. The results show that the BCE method can maintain a second-order accuracy by eliminating the staircase approximation error effectively, whereas the SA and AM methods can only achieve a first-order accuracy. The results of the BCE method show better accuracy than that of the SPS method by considering the physical location of curved surfaces according to the surface location and intersection ratio. In addition, the results of the BCE method agree well with those of the analytical solution method for different dielectric contrast cases.

Reference
1Warnick K FChew W C 2001 Waves Random Media 11 R1
2Jin LZhu Q YFu Y QYu W X 2013 Chin. Phys. 22 104101
3Li JGuo L XZeng HHan X B 2009 Chin. Phys. 18 2757
4Li Q BWu R XYang YSun H L 2013 Chin. Phys. Lett. 30 074208
5Liang CXu Y MWang Z J 2008 Chin. Phys. Lett. 25 3712
6Zhan C LRen X FHuang Y FDuan K MGuo G C 2008 Chin. Phys. Lett. 25 0559
7Lou STsang LChan CIshimaru A 1990 Microw. Opt. Technol. Lett. 3 150
8Hastings F DSchneider J BBroschat S L 1995 IEEE Trans. Antennas Propagal. 43 1183
9Lou STsang LChan C 1991 Waves Random Media 1 287
10Krause KLou STsang LChan C 1991 Microw. Opt. Technol. Lett. 4 255
11Devayya RWingham D 1992 IEEE Trans. Geosci. Remote Sensing 30 645
12Tsang LLou SChan C 1991 Microw. Opt. Technol. Lett. 4 527
13Taflove AHagness S C2005Computational Electrodynamics(Artech House)
14Elsherbeni A ZDemir V2009The Finite-difference Time-domain Method for Electromagnetics with MATLAB SimulationsRaleigh, NCSciTech Pub
15Liu S BLiu S Q 2004 Chin. Phys. 13 1892
16Shen JSha WHuang Z XChen M SWu X L2012Acta Phys. Sin.61190202(in Chinese)
17Liu Y WChen Y WZhang PLiu Z X 2014 Chin. Phys. 23 124102
18Schneider J BShlager K L 1997 IEEE Trans. Antennas Propagal. 45 1830
19Mohammadi ANadgaran HAgio M2005Opt. Express13103671
20Hao YRailton C J 1998 IEEE Trans. Microwave Theory Tech. 46 82
21Fusco M 1990 IEEE Trans. Antennas Propagal. 38 76
22Jurgens TTaflove AUmashankar KMoore T 1992 IEEE Trans. Antennas Propagal. 40 357
23Railton C JCraddock ISchneider J B 1995 Electron. Lett. 31 1585
24Kosmanis TTsiboukis T D 2003 IEEE Trans. Microwave Theory Tech. 51 839
25Oskooi A FKottke CJohnson S G 2009 Opt. Lett. 34 2778
26Liu D CChang H C 2012 IEEE Trans. Antennas Propagal. 60 5259
27Farjadpour ARoundy DRodriguez AIbanescu MBermel PJoannopoulos JJohnson S GBurr G 2006 Opt. Lett. 31 2972
28Kottke CFarjadpour AJohnson S G 2008 Phys. Rev. 77 036611
29Ditkowski ADridi KHesthaven J S 2001 J. Comput. Phys. 170 39
30Xiao TLiu Q H 2004 IEEE Trans. Antennas Propagal. 52 730
31Jackson J D1962Classical ElectrodynamicsNew YorkWiley
32Hwang K PCangellaris A C 2001 IEEE Microw. Wirel. Compon. Lett. 11 158
33Balanis C A2012Advanced Engineering ElectromagneticsVol. 111Wiley Online Library
34Bourlier CPinel NKubické G2013Method of Moments for 2D Scattering Problems: Basic Concepts and ApplicationsJohn Wiley & Sons
35Harrington R F1961Time-harmonic Electromagnetic FieldsMcGraw-Hill
36Kerker M1969The Scattering of LightNew YorkAcademic Press
37Bohren C FHuffman D R2008Absorption and Scattering of Light by Small ParticlesJohn Wiley & Sons